library(tidyverse)
library(haven)
library(AER)
library(Amelia)
library(MatchIt)
library(Matching)
library(MatchingFrontier)
library(ebal)
library(xtable)

#first-stage with geographic filter
####CORRESPONDS TO FIGURE 4a IN PAPER

setwd("C:/Users/dapon/Dropbox/Gov2001 Replication/Original replication download/usedata")
s_filter <- read_dta("survey_data.dta")
survey_data <- read_dta("survey_data.dta")
median <- quantile(unique(survey_data$logdist), na.rm = T)[3] #median(unique(s_filter$logdist),na.rm = T)
s_filter <- s_filter %>% filter(logdist < median)
s_filter <- s_filter %>% filter(!is.na(id) ) 
s_filter <- s_filter %>% filter( !is.na(munid) ) %>% group_by(munid) %>% 
  summarize(treatarrivals2=unique(treatarrivals2), distance=unique(distance), 
            count=n(), treatment=unique(treatment) )

labs <- filter(s_filter, count > 100) %>% 
  mutate(munid_name="",
         munid_name=ifelse(munid==43, "Kos", munid_name), 
         munid_name=ifelse(munid==47, "Lesvos", munid_name), 
         munid_name=ifelse(munid==70, "Samos", munid_name),
         munid_name=ifelse(munid==91, "Chios", munid_name) )


########################
# BINARY TREATMENT MAP #
########################

theme_base <- 
  theme_minimal(base_size=9)  + 
  theme(legend.position=c(0.8,0.8), 
        legend.key.size = unit(0.5, "cm"), 
        axis.text=element_text(size=10))

plt_filter <- ggplot(data=s_filter, aes(x=distance, y=treatarrivals2)) + 
  scale_x_log10()  + theme_base + 
  geom_smooth(aes(col="LOESS"), se=TRUE, method='loess', linetype="dashed") + 
  geom_smooth(aes(col="OLS"), se=TRUE, method='lm', linetype="solid") + 
  geom_point(aes(size=count)) + ylab("Refugee arrivals per capita") + 
  xlab("Distance to Turkey (in km)") + 
  ggtitle("First-Stage Results for Distance < Median") + 
  scale_size(guide=guide_legend(title="Sample size"), breaks=c(50,100,200,300)) + 
  scale_colour_manual(name="Smoothing", values=c("#d7191c", "#2b83ba"), 
                      guide=guide_legend(override.aes=list(fill=NA, cex=0.5, linetype=c("dashed","solid")))) 

plt_filter

setwd("C:/Users/dapon/Dropbox/Gov2001 Replication/Original replication download/extension")

ggsave("firststage_filter_median.pdf",  plt_filter, width=(6.3228344444445/3)*2.5, 
       height=(6.3228344444445/3)*2.5, units="in")


#second-stage with geographic filter
#CORRESPONDS TO FIGURE 4b IN PAPER
setwd("C:/Users/dapon/Dropbox/Gov2001 Replication/Original replication download/usedata")
survey_data <- read_dta("survey_data.dta")
survey_data$mi_m <- survey_data$`_mi_m`
median <- median(unique(survey_data$logdist),na.rm = T)
reg_survey_data <- survey_data %>% filter(logdist < median)

holder_template <- data.frame(score_asylum = rep(NA, 4), # fewer asylum seekers
                              asylumspec_kids = rep(NA,4), # ban from schools 
                              asylumgen = rep(NA,4),  # are a burden
                              asylumspec_terrorism = rep(NA,4), # more crimes
                              asylumspec_criminality  = rep(NA,4),  # more terror attacs
                              asylumspec_burden = rep(NA,4), # component 
                              score_immig = rep(NA, 4), # immigration component
                              incrdecr_economicimmigr = rep(NA, 4),  # fewer economic migrants
                              borders_bordercontrol = rep(NA,4),  # increase border protection
                              diff_muslim = rep(NA,4),  # T(Muslim - Muslim immi)
                              diff_x = rep(NA,4), # T(Christ - Christ immi)
                              score_muslim = rep(NA, 4), # Muslim component
                              incrdecr_minority = rep(NA,4), # decrease representation
                              diff_immig = rep(NA, 4), # T (Christ immi - Muslim immi)
                              islam_acquainted = rep(NA, 4), # how many do not integrate
                              islam_extremism = rep(NA, 4), # how many support extremists
                              score_behavioral = rep(NA, 4), #Behavioral component
                              notifyMPs_rec = rep(NA, 4), # Notify MP? 
                              donation_bin = rep(NA, 4), # Donate?
                              donation_self = rep(NA, 4), # 100-donation 
                              petition_rec = rep(NA, 4), # Sign Petition?
                              score_all = rep(NA, 4), # Overall component
                              diff_nat = rep(NA, 4)# 
)

holder_figure_geo <- holder_template

varlist <- c(  "score_asylum",
               "asylumspec_kids" ,
               "asylumgen" ,
               "asylumspec_terrorism",
               "asylumspec_criminality" ,
               "asylumspec_burden" ,
               "score_immig", 
               "incrdecr_economicimmigr" ,
               "borders_bordercontrol" ,
               "diff_muslim", 
               "diff_x" ,
               "score_muslim" ,
               "incrdecr_minority" ,
               "diff_immig" ,
               "islam_acquainted", 
               "islam_extremism" ,
               "score_behavioral" ,
               "notifyMPs_rec",
               "donation_bin" ,
               "donation_self" ,
               "petition_rec",
               "score_all", 
               "diff_nat")

for ( j in varlist ) { 
  b.out<-NULL
  se.out<-NULL
  for(i in 1:5) {
    reg.data <- reg_survey_data %>% filter(mi_m == i) 
    iv.out <- ivreg(data = reg.data, reg.data[[j]] ~ treatment + female + age + as.factor(edu)| 
                      logdist + female + age + as.factor(edu), weights = weight2_trim)
    cluster.se <- sqrt(diag(vcovCL(iv.out, cluster = ~munid, type = "HC0",
                                   target = NULL, inverse_var = FALSE)))
    b.out <- rbind(b.out, iv.out$coef)
    se.out <- rbind(se.out, cluster.se)
  }
  combined.results <- mi.meld(q = b.out, se = se.out)
  combined.results
  
  holder_figure_geo[1,j] <- combined.results$q.mi[2]
  holder_figure_geo[2,j] <- combined.results$se.mi[2]
  holder_figure_geo[3,j] <- combined.results$q.mi[1]
  holder_figure_geo[4,j] <- combined.results$se.mi[1]
}

fig.geo.filter.data <- as.data.frame(t(holder_figure_geo))
colnames(fig.geo.filter.data) <- c("beta_1", "se_1", "beta_0", "se_0")

fig.geo.filter.data$lo <- fig.geo.filter.data$beta_1 + fig.geo.filter.data$se_1 * -1.96 
fig.geo.filter.data$hi <- fig.geo.filter.data$beta_1 + fig.geo.filter.data$se_1 * 1.96 

fig.geo.filter.data$vars <- rev(c( "Asylum-seeker component (SD=1)",
                                   "Ban from schools (1-5)" ,
                                   "Fewer asylum-seekers (1-5)" ,
                                   "More terror attacks (1-5)",
                                   "More crimes (1-5)" ,
                                   "Are a burden (1-5)" ,
                                   "Immigrant component (SD=1)", 
                                   "Fewer economic migrants (1-5)" ,
                                   "Increase border protection (1-5)" ,
                                   "T(Muslim-Muslim immi.) (1-5)", 
                                   "T(Christ.-Christ. immi.) (1-5)" ,
                                   "Muslim Component (SD=1)" ,
                                   "Decrease representation (1-5)" ,
                                   "T(Christ. immi.-Muslim immi.) (1-5)" ,
                                   "How many do not integrate (1-5)",
                                   "How many support extremists (1-5)" ,
                                   "Behavioral component (SD=1)" ,
                                   "Notify MP? (-2,2)",
                                   "Donate? (0,1)" ,
                                   "(100-donation)/100 (0-1)" ,
                                   "Sign petition? (0,1)",
                                   "Total Component", 
                                   "T(Christ - Muslim) (1-5)"))
fig.geo.filter.data$color <- c(1,2,2,3,3,3,1,2,2,3,3,1,2,3,3,3,1,2,2,2,2,1,1) # gives the right colors

fig.geo.filter.data$vars <- as.character(fig.geo.filter.data$vars) #These lines help order the var list for ggplot
fig.geo.filter.data$vars <- rev(factor(fig.geo.filter.data$vars, levels=unique(fig.geo.filter.data$vars)))

setwd("C:/Users/dapon/Dropbox/Gov2001 Replication/Original replication download/extension")

fig.geo.filter.plot <- ggplot(data = fig.geo.filter.data, aes(x = vars, y = beta_1, color = as.factor(color))) + 
  geom_point() +
  theme_minimal()+ 
  geom_linerange(aes(ymin = lo, ymax = hi)) + coord_flip() + 
  theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.5)) + 
  theme(legend.position="none")+
  theme(plot.title = element_text(size = 8)) + 
  geom_hline(yintercept = 0, linetype = "dotted")+
  scale_color_manual(breaks = c("1", "2", "3"), 
                     values = c("#000000", "#5599ff","#00bb00"))+
  ylab("Logdist < Median, IV")+
  xlab("")  


ggsave(plot=fig.geo.filter.plot, "IV_filter_median.pdf",
       width=(6.3228344444445), height=(6.3228344444445/4)*3, units="in")


#FIRST-STAGE WITH 75TH PERCENTILE CUTOFF
#CORRESPONDS TO FIGURE 5a IN PAPER
setwd("C:/Users/dapon/Dropbox/Gov2001 Replication/Original replication download/usedata")
s_filter <- read_dta("survey_data.dta")
survey_data <- read_dta("survey_data.dta")
third <- quantile(unique(survey_data$logdist), na.rm = T)[4] #median(unique(s_filter$logdist),na.rm = T)
s_filter <- s_filter %>% filter(logdist < third)
s_filter <- s_filter %>% filter(!is.na(id) ) 
s_filter <- s_filter %>% filter( !is.na(munid) ) %>% group_by(munid) %>% 
  summarize(treatarrivals2=unique(treatarrivals2), distance=unique(distance), 
            count=n(), treatment=unique(treatment) )

labs <- filter(s_filter, count > 100) %>% 
  mutate(munid_name="",
         munid_name=ifelse(munid==43, "Kos", munid_name), 
         munid_name=ifelse(munid==47, "Lesvos", munid_name), 
         munid_name=ifelse(munid==70, "Samos", munid_name),
         munid_name=ifelse(munid==91, "Chios", munid_name) )


########################
# BINARY TREATMENT MAP #
########################

theme_base <- 
  theme_minimal(base_size=9)  + 
  theme(legend.position=c(0.8,0.8), 
        legend.key.size = unit(0.5, "cm"), 
        axis.text=element_text(size=10))

plt_filter_third <- ggplot(data=s_filter, aes(x=distance, y=treatarrivals2)) + 
  scale_x_log10()  + theme_base + 
  geom_smooth(aes(col="LOESS"), se=TRUE, method='loess', linetype="dashed") + 
  geom_smooth(aes(col="OLS"), se=TRUE, method='lm', linetype="solid") + 
  geom_point(aes(size=count)) + ylab("Refugee arrivals per capita") + 
  xlab("Distance to Turkey (in km)") + 
  ggtitle("First-Stage Results for Logdist < 75th Percentile") + 
  scale_size(guide=guide_legend(title="Sample size"), breaks=c(50,100,200,300)) + 
  scale_colour_manual(name="Smoothing", values=c("#d7191c", "#2b83ba"), 
                      guide=guide_legend(override.aes=list(fill=NA, cex=0.5, linetype=c("dashed","solid")))) 

plt_filter_third

setwd("C:/Users/dapon/Dropbox/Gov2001 Replication/Original replication download/extension")

ggsave("firstage_filter_75th.pdf",  plt_filter_third, width=(6.3228344444445/3)*2.5, 
       height=(6.3228344444445/3)*2.5, units="in")


#SECOND-STAGE WITH 75TH PERCENTILE CUTOFF
#CORRESPONDS TO FIGURE 5b IN PAPER
setwd("C:/Users/dapon/Dropbox/Gov2001 Replication/Original replication download/usedata")
survey_data <- read_dta("survey_data.dta")
survey_data$mi_m <- survey_data$`_mi_m`
third <- quantile(unique(survey_data$logdist), na.rm = T)[4] #median(unique(s_filter$logdist),na.rm = T)
reg_survey_data <- survey_data %>% filter(logdist < third)

holder_template <- data.frame(score_asylum = rep(NA, 4), # fewer asylum seekers
                              asylumspec_kids = rep(NA,4), # ban from schools 
                              asylumgen = rep(NA,4),  # are a burden
                              asylumspec_terrorism = rep(NA,4), # more crimes
                              asylumspec_criminality  = rep(NA,4),  # more terror attacs
                              asylumspec_burden = rep(NA,4), # component 
                              score_immig = rep(NA, 4), # immigration component
                              incrdecr_economicimmigr = rep(NA, 4),  # fewer economic migrants
                              borders_bordercontrol = rep(NA,4),  # increase border protection
                              diff_muslim = rep(NA,4),  # T(Muslim - Muslim immi)
                              diff_x = rep(NA,4), # T(Christ - Christ immi)
                              score_muslim = rep(NA, 4), # Muslim component
                              incrdecr_minority = rep(NA,4), # decrease representation
                              diff_immig = rep(NA, 4), # T (Christ immi - Muslim immi)
                              islam_acquainted = rep(NA, 4), # how many do not integrate
                              islam_extremism = rep(NA, 4), # how many support extremists
                              score_behavioral = rep(NA, 4), #Behavioral component
                              notifyMPs_rec = rep(NA, 4), # Notify MP? 
                              donation_bin = rep(NA, 4), # Donate?
                              donation_self = rep(NA, 4), # 100-donation 
                              petition_rec = rep(NA, 4), # Sign Petition?
                              score_all = rep(NA, 4), # Overall component
                              diff_nat = rep(NA, 4)# 
)

holder_figure_geo <- holder_template

varlist <- c(  "score_asylum",
               "asylumspec_kids" ,
               "asylumgen" ,
               "asylumspec_terrorism",
               "asylumspec_criminality" ,
               "asylumspec_burden" ,
               "score_immig", 
               "incrdecr_economicimmigr" ,
               "borders_bordercontrol" ,
               "diff_muslim", 
               "diff_x" ,
               "score_muslim" ,
               "incrdecr_minority" ,
               "diff_immig" ,
               "islam_acquainted", 
               "islam_extremism" ,
               "score_behavioral" ,
               "notifyMPs_rec",
               "donation_bin" ,
               "donation_self" ,
               "petition_rec",
               "score_all", 
               "diff_nat")

for ( j in varlist ) { 
  b.out<-NULL
  se.out<-NULL
  for(i in 1:5) {
    reg.data <- reg_survey_data %>% filter(mi_m == i) 
    iv.out <- ivreg(data = reg.data, reg.data[[j]] ~ treatment + female + age + as.factor(edu)| 
                      logdist + female + age + as.factor(edu), weights = weight2_trim)
    cluster.se <- sqrt(diag(vcovCL(iv.out, cluster = ~munid, type = "HC0",
                                   target = NULL, inverse_var = FALSE)))
    b.out <- rbind(b.out, iv.out$coef)
    se.out <- rbind(se.out, cluster.se)
  }
  combined.results <- mi.meld(q = b.out, se = se.out)
  combined.results
  
  holder_figure_geo[1,j] <- combined.results$q.mi[2]
  holder_figure_geo[2,j] <- combined.results$se.mi[2]
  holder_figure_geo[3,j] <- combined.results$q.mi[1]
  holder_figure_geo[4,j] <- combined.results$se.mi[1]
}

fig.geo.filter.data <- as.data.frame(t(holder_figure_geo))
colnames(fig.geo.filter.data) <- c("beta_1", "se_1", "beta_0", "se_0")

fig.geo.filter.data$lo <- fig.geo.filter.data$beta_1 + fig.geo.filter.data$se_1 * -1.96 
fig.geo.filter.data$hi <- fig.geo.filter.data$beta_1 + fig.geo.filter.data$se_1 * 1.96 

fig.geo.filter.data$vars <- rev(c( "Asylum-seeker component (SD=1)",
                                   "Ban from schools (1-5)" ,
                                   "Fewer asylum-seekers (1-5)" ,
                                   "More terror attacks (1-5)",
                                   "More crimes (1-5)" ,
                                   "Are a burden (1-5)" ,
                                   "Immigrant component (SD=1)", 
                                   "Fewer economic migrants (1-5)" ,
                                   "Increase border protection (1-5)" ,
                                   "T(Muslim-Muslim immi.) (1-5)", 
                                   "T(Christ.-Christ. immi.) (1-5)" ,
                                   "Muslim Component (SD=1)" ,
                                   "Decrease representation (1-5)" ,
                                   "T(Christ. immi.-Muslim immi.) (1-5)" ,
                                   "How many do not integrate (1-5)",
                                   "How many support extremists (1-5)" ,
                                   "Behavioral component (SD=1)" ,
                                   "Notify MP? (-2,2)",
                                   "Donate? (0,1)" ,
                                   "(100-donation)/100 (0-1)" ,
                                   "Sign petition? (0,1)",
                                   "Total Component", 
                                   "T(Christ - Muslim) (1-5)"))
fig.geo.filter.data$color <- c(1,2,2,3,3,3,1,2,2,3,3,1,2,3,3,3,1,2,2,2,2,1,1) # gives the right colors

fig.geo.filter.data$vars <- as.character(fig.geo.filter.data$vars) #These lines help order the var list for ggplot
fig.geo.filter.data$vars <- rev(factor(fig.geo.filter.data$vars, levels=unique(fig.geo.filter.data$vars)))

setwd("C:/Users/dapon/Dropbox/Gov2001 Replication/Original replication download/extension")

fig.geo.filter.plot <- ggplot(data = fig.geo.filter.data, aes(x = vars, y = beta_1, color = as.factor(color))) + 
  geom_point() +
  theme_minimal()+ 
  geom_linerange(aes(ymin = lo, ymax = hi)) + coord_flip() + 
  theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.5)) + 
  theme(legend.position="none")+
  theme(plot.title = element_text(size = 8)) + 
  geom_hline(yintercept = 0, linetype = "dotted")+
  scale_color_manual(breaks = c("1", "2", "3"), 
                     values = c("#000000", "#5599ff","#00bb00"))+
  ylab("Logdist < 75th Percentile, IV")+
  xlab("")  

ggsave(plot=fig.geo.filter.plot, "IV_filter_75.pdf", 
       width=(6.3228344444445), height=(6.3228344444445/4)*3, units="in")




#FIRST-STAGE WITH 75TH PERCENTILE CUTOFF
#CORRESPONDS TO FIGURE 6a IN PAPER
setwd("C:/Users/dapon/Dropbox/Gov2001 Replication/Original replication download/usedata")
s_filter <- read_dta("survey_data.dta")
survey_data <- read_dta("survey_data.dta")
first <- quantile(unique(survey_data$logdist), na.rm = T)[2] #median(unique(s_filter$logdist),na.rm = T)
s_filter <- s_filter %>% filter(logdist < first)
s_filter <- s_filter %>% filter(!is.na(id) ) 
s_filter <- s_filter %>% filter( !is.na(munid) ) %>% group_by(munid) %>% 
  summarize(treatarrivals2=unique(treatarrivals2), distance=unique(distance), 
            count=n(), treatment=unique(treatment) )

labs <- filter(s_filter, count > 100) %>% 
  mutate(munid_name="",
         munid_name=ifelse(munid==43, "Kos", munid_name), 
         munid_name=ifelse(munid==47, "Lesvos", munid_name), 
         munid_name=ifelse(munid==70, "Samos", munid_name),
         munid_name=ifelse(munid==91, "Chios", munid_name) )


########################
# BINARY TREATMENT MAP #
########################

theme_base <- 
  theme_minimal(base_size=9)  + 
  theme(legend.position=c(0.8,0.8), 
        legend.key.size = unit(0.5, "cm"), 
        axis.text=element_text(size=10))

plt_filter_first <- ggplot(data=s_filter, aes(x=distance, y=treatarrivals2)) + 
  scale_x_log10()  + theme_base + 
  geom_smooth(aes(col="LOESS"), se=TRUE, method='loess', linetype="dashed") + 
  geom_smooth(aes(col="OLS"), se=TRUE, method='lm', linetype="solid") + 
  geom_point(aes(size=count)) + ylab("Refugee arrivals per capita") + 
  xlab("Distance to Turkey (in km)") + 
  ggtitle("First-Stage Results for Logdist < 25th Percentile") + 
  scale_size(guide=guide_legend(title="Sample size"), breaks=c(50,100,200,300)) + 
  scale_colour_manual(name="Smoothing", values=c("#d7191c", "#2b83ba"), 
                      guide=guide_legend(override.aes=list(fill=NA, cex=0.5, linetype=c("dashed","solid")))) 

plt_filter_first

setwd("C:/Users/dapon/Dropbox/Gov2001 Replication/Original replication download/extension")

ggsave("firstage_filter_25th.pdf",  plt_filter_first, width=(6.3228344444445/3)*2.5, 
       height=(6.3228344444445/3)*2.5, units="in")


#SECOND-STAGE WITH 25TH PERCENTILE CUTOFF
#CORRESPONDS TO FIGURE 6b IN PAPER
setwd("C:/Users/dapon/Dropbox/Gov2001 Replication/Original replication download/usedata")
survey_data <- read_dta("survey_data.dta")
survey_data$mi_m <- survey_data$`_mi_m`
first <- quantile(unique(survey_data$logdist), na.rm = T)[2] #median(unique(s_filter$logdist),na.rm = T)
reg_survey_data <- survey_data %>% filter(logdist < first)

holder_template <- data.frame(score_asylum = rep(NA, 4), # fewer asylum seekers
                              asylumspec_kids = rep(NA,4), # ban from schools 
                              asylumgen = rep(NA,4),  # are a burden
                              asylumspec_terrorism = rep(NA,4), # more crimes
                              asylumspec_criminality  = rep(NA,4),  # more terror attacs
                              asylumspec_burden = rep(NA,4), # component 
                              score_immig = rep(NA, 4), # immigration component
                              incrdecr_economicimmigr = rep(NA, 4),  # fewer economic migrants
                              borders_bordercontrol = rep(NA,4),  # increase border protection
                              diff_muslim = rep(NA,4),  # T(Muslim - Muslim immi)
                              diff_x = rep(NA,4), # T(Christ - Christ immi)
                              score_muslim = rep(NA, 4), # Muslim component
                              incrdecr_minority = rep(NA,4), # decrease representation
                              diff_immig = rep(NA, 4), # T (Christ immi - Muslim immi)
                              islam_acquainted = rep(NA, 4), # how many do not integrate
                              islam_extremism = rep(NA, 4), # how many support extremists
                              score_behavioral = rep(NA, 4), #Behavioral component
                              notifyMPs_rec = rep(NA, 4), # Notify MP? 
                              donation_bin = rep(NA, 4), # Donate?
                              donation_self = rep(NA, 4), # 100-donation 
                              petition_rec = rep(NA, 4), # Sign Petition?
                              score_all = rep(NA, 4), # Overall component
                              diff_nat = rep(NA, 4)# 
)

holder_figure_geo <- holder_template

varlist <- c(  "score_asylum",
               "asylumspec_kids" ,
               "asylumgen" ,
               "asylumspec_terrorism",
               "asylumspec_criminality" ,
               "asylumspec_burden" ,
               "score_immig", 
               "incrdecr_economicimmigr" ,
               "borders_bordercontrol" ,
               "diff_muslim", 
               "diff_x" ,
               "score_muslim" ,
               "incrdecr_minority" ,
               "diff_immig" ,
               "islam_acquainted", 
               "islam_extremism" ,
               "score_behavioral" ,
               "notifyMPs_rec",
               "donation_bin" ,
               "donation_self" ,
               "petition_rec",
               "score_all", 
               "diff_nat")

for ( j in varlist ) { 
  b.out<-NULL
  se.out<-NULL
  for(i in 1:5) {
    reg.data <- reg_survey_data %>% filter(mi_m == i) 
    iv.out <- ivreg(data = reg.data, reg.data[[j]] ~ treatment + female + age + as.factor(edu)| 
                      logdist + female + age + as.factor(edu), weights = weight2_trim)
    cluster.se <- sqrt(diag(vcovCL(iv.out, cluster = ~munid, type = "HC0",
                                   target = NULL, inverse_var = FALSE)))
    b.out <- rbind(b.out, iv.out$coef)
    se.out <- rbind(se.out, cluster.se)
  }
  combined.results <- mi.meld(q = b.out, se = se.out)
  combined.results
  
  holder_figure_geo[1,j] <- combined.results$q.mi[2]
  holder_figure_geo[2,j] <- combined.results$se.mi[2]
  holder_figure_geo[3,j] <- combined.results$q.mi[1]
  holder_figure_geo[4,j] <- combined.results$se.mi[1]
}

fig.geo.filter.data <- as.data.frame(t(holder_figure_geo))
colnames(fig.geo.filter.data) <- c("beta_1", "se_1", "beta_0", "se_0")

fig.geo.filter.data$lo <- fig.geo.filter.data$beta_1 + fig.geo.filter.data$se_1 * -1.96 
fig.geo.filter.data$hi <- fig.geo.filter.data$beta_1 + fig.geo.filter.data$se_1 * 1.96 

fig.geo.filter.data$vars <- rev(c( "Asylum-seeker component (SD=1)",
                                   "Ban from schools (1-5)" ,
                                   "Fewer asylum-seekers (1-5)" ,
                                   "More terror attacks (1-5)",
                                   "More crimes (1-5)" ,
                                   "Are a burden (1-5)" ,
                                   "Immigrant component (SD=1)", 
                                   "Fewer economic migrants (1-5)" ,
                                   "Increase border protection (1-5)" ,
                                   "T(Muslim-Muslim immi.) (1-5)", 
                                   "T(Christ.-Christ. immi.) (1-5)" ,
                                   "Muslim Component (SD=1)" ,
                                   "Decrease representation (1-5)" ,
                                   "T(Christ. immi.-Muslim immi.) (1-5)" ,
                                   "How many do not integrate (1-5)",
                                   "How many support extremists (1-5)" ,
                                   "Behavioral component (SD=1)" ,
                                   "Notify MP? (-2,2)",
                                   "Donate? (0,1)" ,
                                   "(100-donation)/100 (0-1)" ,
                                   "Sign petition? (0,1)",
                                   "Total Component", 
                                   "T(Christ - Muslim) (1-5)"))
fig.geo.filter.data$color <- c(1,2,2,3,3,3,1,2,2,3,3,1,2,3,3,3,1,2,2,2,2,1,1) # gives the right colors

fig.geo.filter.data$vars <- as.character(fig.geo.filter.data$vars) #These lines help order the var list for ggplot
fig.geo.filter.data$vars <- rev(factor(fig.geo.filter.data$vars, levels=unique(fig.geo.filter.data$vars)))

setwd("C:/Users/dapon/Dropbox/Gov2001 Replication/Original replication download/extension")

fig.geo.filter.plot <- ggplot(data = fig.geo.filter.data, aes(x = vars, y = beta_1, color = as.factor(color))) + 
  geom_point() +
  theme_minimal()+ 
  geom_linerange(aes(ymin = lo, ymax = hi)) + coord_flip() + 
  theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.5)) + 
  theme(legend.position="none")+
  theme(plot.title = element_text(size = 8)) + 
  geom_hline(yintercept = 0, linetype = "dotted")+
  scale_color_manual(breaks = c("1", "2", "3"), 
                     values = c("#000000", "#5599ff","#00bb00"))+
  ylab("Logdist < 25th Percentile, IV")+
  xlab("")  

ggsave(plot=fig.geo.filter.plot, "IV_filter_25.pdf", width=(6.3228344444445), height=(6.3228344444445/4)*3, units="in")

